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NUMERICAL  METHODS  FOR  SINGULAR  PERTURBATION  PROBLEMS 


INTRODUCTION 


Numerical  analysis  of  singular  perturbation  problems  has  been  of 
considerable  interest  recently  since  these  problems  mimic  difficulties 
encountered  with  the  Navier-Stokes  equations.  The  existence  of  a  small 
parameter,  e,  in  singular  perturbation  problems  causes  the  appearance  of  a 
thin  boundary  layer  in  their  solutions.  If  £  tends  to  zero,  an  infinite 
number  of  grid  points  need  to  be  clustered  near  the  physical  boundary  so 
that  the  boundary  layer  can  be  properly  resolved.  Moreover,  in  the 
nonlinear  case,  the  location  of  the  boundary  layer  is  not  known  a  priori, 
and  this  renders  the  problem  more  difficult. 

Finite  difference  approximations  of  singular  perturbation  problems  lead 
to  unsynmetric  matrices.  Reliable  and  fast  iterative  methods  will  be 
presented  here  for  linear  and  nonlinear  problems.  The  assumption  is  that 
matrices  are  diagonally  dominant,  and  the  methods  do  not  require  the 
symmetry  property  of  matrices.  Numerical  results  will  be  compared  with 
solutions  from  standard  algorithms. 


THE  LINEAR  PROBLEM 


Consider  the  linear  advection-diffusion  equation  in  conservative  form: 


where  0  is  the  advected  quantity,  u  and  v  are  velocity  components  in  x  and 
y  directions  respectively,  and  e  is  the  diffusivity  constant.  Fiadeiro  and 


1 


Veronis^  proposed  the  weighted  mean  (W-M)  scheme  to  discretize  equation 
(1)  as 

M  <j) .  .  .  +  E  <J>.  -  .  +  S  ,  +  N  ,  -  C  <j>.  .  *  0  , 

i“l,J  *1+1, J  i.J-1  i,J+l  i,J  ’ 
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s  =  J.  -.1/2 
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N  2Ay 


coth 


coth 


coth 


coth 


-  1/2  Ax 
2e 


1/2  Ax 


+  1 


C=W+E+S+N. 


The  W-M  scheme  retains  the  five-point  operator  character  of  central 
difference  except  that,  in  the  limit,  it  tends  to  central  difference  for 
strongly  diffusive  cases  and  to  upwind  difference  for  strongly  advective 
cases.  This  scheme  is  derived  by  using  the  flux  conserving  principle;  the 
discretization  error  is  0  (Ax^,  Ay^).  Significant  advantages  of  the 
scheme  are  that  there  are  no  local  instabilities  in  the  solution  with 
respect  to  central  difference  and  that  no  testing  is  required  at  each  grid 
point  with  respect  to  upwind  difference.  Its  stability  for  all  grid 
spacings  has  been  found  to  result  in  fast  convergence  when  using  the 
preconditioned  minimal  residual  method  and  the  multigrid  method.  For 
comparative  purposes,  the  successive  line  over-relaxation  (SLOR)  algorithm 
which  serves  as  a  standard  iterative  method,  is  briefly  described. 
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SUCCESSIVE  LINE  OVER- RELAXATION  (SLOR)  METHOD 


Equation  (1)  can  be  rewritten  as 
L<p  -  0, 

2 

where  L  is  a  differential  operator.  The  SLOR  method  can  be  implemented 
as  follows: 


a  +  6  + 

yy 


n 
e .  . 


-  0!  L 


1J 


where  e?.  =  -  d>?.  and  cj  is  the  relaxation  parameter. 

The  scheme  is  implicit  in  the  sense  that  in  the  x  and  y  directions  an 
inversion  of  a  tridiagonal  matrix  equation  for  a  given  value  of  i  or  j  is 
required.  In  these  numerical  examples  an  optimal  relaxation  factor,  U),  is 
assumed. 


THE  PRECONDITIONED  MINIMAL  RESIDUAL  (PMR)  METHOD 

The  discretization  of  equation  (1)  yields  a  large  system  of  linear 
equations 

Ax  *  b,  (2) 

where  A  is  a  large,  sparse,  and  nonsymmetric  coefficient  matrix;  x  is  the 
solution  vector;  and  b  contains  boundary  terms. 

Introduced  now  is  a  nonsingular  matrix,  C,  called  the  preconditioning 
matrix.  Equation  (2)  is  thus  equivalent  to 

AA  A 

Ax  *  b,  (3) 

^  ^  /A 

where  A  *  C_1A,  x  ■  x  and  b  ■  C-1b.  This  case  is  denoted  as  PMR  I.  A  more 
general  choice  for  equation  (2)  is  A  m  PAQ,  X  *  Q~^x,  and  b  ■  Pb,  where 
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P  and  Q  are  nonsingular  matrices.  The  second  case  is  called  PMR  II.  The 

procedure  for  transforming  equation  (2)  into  equation  (3)  is  known  as  a 

preconditioning  technique.  Once  the  matrix  C  or  the  matrices  P  and  Q  have 

3 

been  chosen,  the  standard  minimal  residual  (MR)  method  can  be  applied  to 
the  transformed  equation.  A  different  choice  for  C  or  for  P  and  Q,  there¬ 
fore,  leads  to  a  different  preconditioned  MR  method  and,  in  consequence, 
different  rates  of  convergence. 

To  accelerate  the  convergence  rate  of  the  PMR  method,  it  is  desirable  to 
find  matrices  such  that  K(C  ^A)  «  K(A)  or  K(PAQ)  «  K(A),  where  K(M)  is 
the  condition  number  of  matrix  M.  The  PMR  I  method  can  be  implemented  by 
the  following  algorithm: 

Let  K  =  0  ,  rQ  »  b  *  A  X.  For  K=0,.l,  2,  ...  ,  compute  the 

A 

vectors  *K+1  f  rom 

\  rK  '  “he”  “k  ■  ’  tK  )  ’ 

rK+l  =  ^  “  A  *K+1  rK  "  °K  ^rK  ’ 

A  similar  procedure  can  be  applied  to  the  PMR  II  method.  To  minimize 
the  computational  work  involved  in  the  preconditioning  technique,  the 
matrix  C  should  easily  be  invertible.  In  general,  C  is  chosen  to  be  the 
product  of  lower  and  upper  triangular  matrices;  i.e.,  C  ■  LU.  The  inverse 
can  therefore  be  obtained  through  simple  forward  and  backward  substitutions. 
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ROW-SUM  AGREEMENT  FACTORIZATION 


Wong  suggested  an  incomplete  factorization  based  on  the  principle  of 
agreement  of  row-sum  between  the  preconditioning  matrix  C  and  the  original 
matrix  A: 

Row-sum  (A)  ■  Row-sum  (C)  ■  Row-sum  (LU) 

A  ■  LU+R,  where  R  is  the  defect  matrix;  hence  Row-sum  (R)  =  0 
Therefore,  matrix  C  resembles  matrix  A  in  the  sense  that  their 
eigenvalues  have  the  same  bound,  namely  the  norm  of  A.  Different  r  irithms 

#  Zf 

can  be  found  m  Wong  for  constructing  L  and  U.  The  row-sum  agreem 
factorization  is  shown  schematically  in  figure  1.  L  and  U  have  not 
elements  in  these  positions  which  correspond  to  the  nonzero  elements  j.n  the 
lower  and  upper  triangular  parts  of  A.  The  off-diagonal  elements  of  C, 
whose  locations  correspond  to  the  nonzero  off-diagonal  elements  of  A,  are 
set  to  those  values.  If  ||r||  is  small,  LU  will  be  a  good  approximation  to 
A;  consequently,  a  fast  convergence  can  be  achieved  by  the  PMR  method. 


THE  MULTIGRID  PROCEDURE 


The  «Cycle  C»  algorithm  with  correction  storage  is  used  here,  as 
presented  by  Brandt.^  The  initial  approximate  solution  is  defined  on  the 
finer  grid,  and  then  relaxed.  When  the  relaxation  is  found  to  be 
inefficient,  a  correction  is  projected  on  a  coarser  grid.  When  the 
relaxation  converges  on  any  grid,  the  values  on  the  finer  grid  are  updated, 
and  relaxation  performed  at  this  finer  level.  The  sequence  of  grids  used  is 
not  defined  a  priori,  but  depends  on  the  efficiency  and  convergence  of  the 
relaxation  just  completed.  Data  are  transferred  from  fine  to  coarse  grids 


by  projection,  and  from  coarse  to  fine  grids  by  linear  interpolation.  The 
relaxation  performed  is  point  Gauss-Seidel .  Since  the  problem  is  linear  and 
the  velocity  field  is  explicitly  given,  the  coefficients  N,  S,  E,  W  are 
explicitly  computed  at  each  level,  using  the  appropriate  location  (x,y)  and 
Ax,  Ay. 

Used  here  for  an  efficiency  criterion  is: 

residual  at  present  relaxation  <  0.6  (residual  at  former 
relaxation);  and  for  a  convergence  criterion: 


residual  on  "coarse"  grid  K\  <  0.3  /residual  on  "fine"  grid  K  +  1 


This  scaling  takes  into  account  the  fact  that  the  method  is  second 
2 

order,  i.e.,  error — h  . 

A  simple  local  analysis  produces  an  interesting  result  for  variable  flow 
fields.  Consider  the  smoothing  rate  for  Gauss-Seidel  relaxation  in  the 
usual  ordering,  with  both  x  and  y  increasing: 

p  -  MAX. 

x  <  E,  ,  n  <  n 
2 

Since  p  is  the  reduction  factor  for  high  frequency  modes,  it  should  be  as 
small  as  possible.  But,  from  the  explicit  expression  for  N,  S,  E,  W,  it  can 
be  seen  that  if  u  or  v  is  negative  and  large  in  magnitude,  E  »  W  and  N  >>  S, 
as  they  should  be  in  an  upstream  differentiation  formulation.  This  produces 
a  smoothing  rate  very  near  to  one.  Thus,  if  u  is  negative,  relaxation 
should  not  be  in  the  direction  of  increasing  x. 


7 


If  Che  flow  field  is  reasonably  simple,  so  ChaC  Che  relaxacion  can  be 
explicidy  performed  following  Che  flow,  i.e.,  x,  y,  increasing  if  only  u, 
v  >  0,  Chen  Che  convergence  is  much  faster. 

IC  wasn'c  necessary  Co  use  chis  meChod  on  Chese  cesC  problems  as  u,  v  >  0 
everywhere,  buC  Chere  are  resulCs  showing  speedup  of  convergence  of  a  Benard 
cell  Cype  problem  solved  by  LusCman,  eC  al.^  In  any  case,  mulcigrid 
iceracion  is  evidenCly  faster  chan  SLOR,  even  when  Che  besC  overrelaxation 
parameter  is  used,  as  Cable  1  shows. 


NUMERICAL  RESULTS  OF  A  LINEAR  EXAMPLE 


Consider  Che  following  model  problem: 


where  u  =*  U(1  -  X/L),  v  =  UY/L  ,  and  P  *  UL/K,  with  boundary  conditions: 

$  *  0  at  left  and  bottom, 

4>  *  1  at  right  and  top. 

P  is  the  Peclet  number  which  measures  the  intensity  of  advection  relative  to 
diffusion.  U,  L,  and  K  (diffusivity  constant)  are  quantities  used  for 
nondimensionalizacion.  Figure  2  depicts  Che  flow  pattern  where  streamlines 
enter  along  the  left  boundary  and  leave  through  the  top  boundary.  Plots  of 
isotherms  are  shown  in  figures  3  and  4  for  P  -  1  and  100,  respectively. 
Observe  the  thin  boundary  layer  developed  near  x  *  L  and  y  “  L  for  strong 


advective  cases.  In  table  1  the  number  of  iterations  for  P  from  10  to  400 


Figure  3.  Isotherms  for  P  =  1  Figure  4.  Isotherms  for  P  =  100 
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Table  1.  Number  of  Iterations  with  Different  Relaxation  Methods 


1 


Peclet  No. 

Relaxation 

Method 

10 

— 

50 

100 

— 

400 

SLOR  in  x  and  y  Directions* 

120 

73 

60 

23 

PMR  I 

45 

48 

37 

22 

PMR  II 

93 

45 

35 

21 

Multigrid,  WUi 

95.6 

77.2 

________ 

64.0 

45.2 

■ 

*One  SLOR  iteration  is  a  sweep  in  the  x  direction,  followed  by  a  sweep 
in  the  y  direction.  It  is  approximately  equal  to  the  work  done  by  one  PMR 
iteration. 

tOne  WU  (work  unit)  is  equivalent  to  one  point  relaxation  sweep  over  a 
61  x  61  grid.  It  is  about  half  the  work  done  in  one  SLOR  relaxation. 


is  compared  for  a  grid  61  x  61.  Convergence  was  reached  when  the  norm  of 
the  residual  vector  was  less  than  10  4 

For  all  four  iteration  methods  considered,  it  was  found  that  the  W-M 
scheme  converges  faster  for  more  convective  flows  since  the  preconditioning 
matrix  C  becomes  closer  to  the  original  matrix  A  as  the  Peclet  number 
becomes  larger.  The  PMR  and  multigrid  methods  need  about  the  same  storage 
and  are  efficient  for  high  as  well  as  low  Peclet  numbers.  Furthermore,  no 
parameter  is  required  for  the  PMR  algorithm.  The  multigrid  program  involves 
some  complexity  since  the  various  sizes  of  the  arrays  represent  the  same 
quantity  and  there  is  a  need  for  interpolations.  On  the  other  hand,  the  PMR 
method  needs  an  efficient  LU  decomposition. 


THE  NONLINEAR  PROBLEM 


Consider  the  following  2-D  Burger's  equation  in  a  unit  square  in  R^ 


! 


(x,  y) 


0  <  x  <  1,  0  <  y  <  1 


*  (4  =  *  (4 


+  e  V  u  ■  0  , 


(4) 


with  boundary  conditions: 

u  ■  0  at  left  and  bottom, 
u  ■  1  at  right  and  top. 

For  standard  upwind  differencing  of  the  convective  terms,  the 
nonlinearity  might  cause  instabilities  and  convergence  to  the  incorrect 
(nonphysical)  solution.  Hence,  the  monotone  difference  scheme  of  Engquist 
and  Osher^  is  used  here  to  approximate  the  nonlinear  convective  terms,  and 
central  difference  is  used  to  discretize  the  Laplacian  operator.  For  the 


11 


w 
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1-D  singular  perturbation  equation,  Osher  has  proved,  via  the  time- 

dependent  equation,  that  convergence  to  a  unique  steady-state  solution 

occurs.  Moreover,  computational  evidence  of  the  1-D  case  in  Osher' s  paper 

indicates  that  the  method  is  very  useful  in  accurately  locating  regions  of 

large  gradient  even  for  complicated  problems.  For  nonuniform  grids, 

9 

Abrahams son  and  Osher  have  proved  general  convergence  results  and  have 
shown  that  the  Engquist-Osher  scheme  reproduces  essential  properties  of  the 
true  solution. 

Next,  the  monotone  difference  scheme  is  described  briefly.  Let 
1  2 

f(u)  =  —  u  .  The  Engquist-Osher  scheme  is  based  on  the  approximation 


where 


_f(u)  x~~k  f-  (ui)  +  A-  f+  (ui) 


\  ui  *  ui+l  -  ,  A  -  u.  -  u.  -  U..J 


and 


f+  (ui)  “  f  [max  (ui»  °)  »  f-  (ui)  “  f  min  (ui>  °) 

For  each  grid  point,  the  discretized  form  of  equation  (4)  reads: 

£  u.  .  +  u.  ,  .  +  u.  .  ,  +  u.  .  ,  -  4  u.  .  I 

I  J+l  i»jJ 

-  h  [f  '  -  2  £-  *  2  f. 


'  £*  -  £.  (»if  j-i) 


0, 


(5) 


(6) 
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where  h  is  the  uniform  grid  size  and 


f  (u) 


f+  (u) 


/ 


2 

0  , 


u  <  0 

u  >_  0 

u  >  0 

u  <  0 


Convergence  to  ^steady-state  solution  will  be  accomplished  via  artificial 
time  stepping  and  Newton's  method. 


METHOD  OF  TIME  STEPPING  AND  NEWTON'S  METHOD 


u£  »  ~  u^ ) y  +  e  (7) 

with  boundary  conditions  and  initial  conditions. 

Discretization  of  equation  (7)  with  respect  to  time  yields 


n+1 

u.  . 


A  . 

At  Lu 


where  Lu11  represents  the  right  hand  side  of  equation  (7)  and  Atn  the 
time  increment.  In  this  example,  the  time  steps,  Atn,  which  are  limited 
by  the  Courant-Friedrich-Levy  (CFL)  condition,  are  chosen  to  be  fixed  for  all 
n,  Atn  *  At°  since  the  finite  difference  equation  satisfies  the  maximum 
principle. 

Instead  of  solving  the  time-dependent  equation  (7),  Newton's  method  can 
be  used  to  iterate  an  initial  guess  of  the  steady-state  solution.  However, 


to  ensure  convergence,  the  starting  guess  of  solution  should  be  close  enough 
to  the  exact  solution. 

Consider  now  the  basic  idea  of  iterative  methods  for  nonlinear  systems 
of  equations.  Any  iterative  scheme  can  be  expressed  as 


M6un  =  -  Lun  ,  (8) 

where  6un  =  u11  +  *  -  un,  Lu11  is  the  residual  vector  at  the  n^  iteration, 
and  M  is  an  iterative  matrix.  It  appears  that  many  choices  of  M  are 
possible.  The  closer  that  M  is  to  representing  L,  the  faster  convergence 
is.  In  general,  it  is  required  that  M  will  be  chosen  in  such  a  way  that: 

1.  M  can  be  computed  easily. 

2.  Only  a  little  storage  is  needed  to  store  M. 

3.  (Su11  can  be  computed  rapidly  from  equation  (8). 

4.  The  iteration  scheme  as  expressed  in  equation  (8)  should  converge. 

If  M  is  selected  to  be  the  Jacobian  of  L,  then  equation  (8)  will  be  a  Newton 

iterative  scheme  and  the  rate  of  convergence  will  be  quadratic. 

til 

In  this  example  the  Jacobian  of  L  is  given,  at  n  step,  by 


M6u 


n 


=  e 


+ 


p  n  \  n  p  n 
Ou  I  -  u  6u 

yy  /  x 


n.  n 
u  ou 

y 


Again  the  convective  and  diffusive  parts  are  approximated  by  the  W-M  scheme 
since  equation  (9)  is  linear  in  Sun.  Then  the  calculations  of  the  cor¬ 
rection  vectors,  Su11,  should  invoke  the  combination  of  the  W-M  scheme  and 
the  preconditioned  minimal  residual  method  for  efficiency  reasons.  Note 
that  in  using  the  PMR  method  to  solve  the  linearized  system  of  equation  (8), 
only  a  few  iterations  are  needed  for  each  outer  iteration,  since  the 
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interest  is  in  the  convergence  for  the  whole  nonlinear  problem.  The  number 
of  inner  iterations  is  fixed  at  six  in  these  numerical  examples. 


A  NONLINEAR  EXAMPLE 


The  combined  Newton  and  PMR  methods  discussed  in  the  previous  section 
result  in  a  very  fast  converging  iterative  method  for  the  solution  of 
nonlinear  systems  of  equations  that  approximate  equation  (4).  In  table  2 
the  number  of  iterations  is  compared  for  the  explicit  scheme  and  Newton's 
method  for  e  =  0.1,  0.05,  and  0.01  with  a  grid  31  x  31.  Note  that  all 
calculations  were  started  with  the  same  uniform  initial  guess  of  u  *  0.5, 
and  convergence  was  reached  when  the  norm  of  the  residual  vector  was  less 

-4 

than  10  .  Table  3  shows  the  total  computing  time  on  a  VAX  11/780 

computer  for  the  two  methods.  There  are  slight  changes  in  the  number  of 
Newton's  iterations  and  drastic  increases  in  the  number  of  time  steps  when  € 
decreases.  Table  4  lists  the  number  of  Newton's  iterations  and  the 
computing  time  with  e  =  0.01  for  different  starting  guesses  of  solutions. 


CONCLUSIONS 


The  PMR  and  multigrid  methods  are  very  efficient  iterative  procedures 
for  large  as  well  as  small  e.  Moreover,  numerical  results  indicate  that 
the  stronger  the  convection  is,  the  faster  the  convergence  is.  Finally,  it 
has  been  shown  that  Newton's  method  combined  with  the  PMR  method  can  be 
applied  very  successfully  to  a  nonlinear  multidimensional  problem. 
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Table  2.  Iteration  Numbers 


e 

0.1 

0.05  ! 

0.01 

Time  Steps 

68 

134 

437 

Newton  Iterations 

8  . 

9 

.  _  _  J 

7 

Table  3.  CPU  Time 


e 

1 

0.1 

0.05 

: 

0.01 

Explicit  Scheme 

19.72 

32.01 

1:26.45 

Newton's  Mechod 

31.99 

1— 

33.46 

29.85 

Table  4.  Runs  with  Different  Initial  Guesses 


£  =  0.01 

i 

II 

III 

Newton  Iterations 

9 

7 

8 

Computing  Time 

-  — 

32.94 

!  29.85 

31.51 

I  *  uniform  initial  guess  of  0 

II  *  uniform  initial  guess  of  0.5 

III  ■  uniform  initial  guess  of  1 
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